R Notebook com GWR

Este notebook R realiza uma análise de regressão ponderada geograficamente (GWR) para modelar a relação entre roubos e uma variável preditora de um arquivo shapefile. A GWR é uma técnica estatística que permite modelar relações espaciais não estacionárias, ou seja, relações que variam geograficamente. Ao contrário da regressão linear global, que assume que a relação entre as variáveis é constante em toda a área de estudo, a GWR permite que os coeficientes de regressão variem espacialmente, capturando a heterogeneidade espacial.

Função para verificar e instalar pacotes

install_if_missing <- function(packages) {
  new_packages <- packages[!packages %in% installed.packages()[, "Package"]]
  if (length(new_packages) > 0) {
    install.packages(new_packages)
  }
}

necessary_packages <- c("sf", "spdep", "tidyverse", "ncf", "tmap", "GWmodel", "dplyr")
install_if_missing(necessary_packages)


suppressPackageStartupMessages({
  library(sf)
  library(spdep)
  library(tidyverse)
  library(ncf)
  library(tmap)
  library(GWmodel)
  library(dplyr)
})

## Warning: pacote 'sf' foi compilado no R versão 4.4.3

## Warning: pacote 'spdep' foi compilado no R versão 4.4.3

## Warning: pacote 'spData' foi compilado no R versão 4.4.3

## Warning: pacote 'tidyverse' foi compilado no R versão 4.4.3

## Warning: pacote 'ggplot2' foi compilado no R versão 4.4.3

## Warning: pacote 'tibble' foi compilado no R versão 4.4.3

## Warning: pacote 'tidyr' foi compilado no R versão 4.4.3

## Warning: pacote 'readr' foi compilado no R versão 4.4.3

## Warning: pacote 'dplyr' foi compilado no R versão 4.4.3

## Warning: pacote 'forcats' foi compilado no R versão 4.4.3

## Warning: pacote 'lubridate' foi compilado no R versão 4.4.3

## Warning: pacote 'tmap' foi compilado no R versão 4.4.3

## Warning: pacote 'GWmodel' foi compilado no R versão 4.4.3

## Warning: pacote 'robustbase' foi compilado no R versão 4.4.3

## Warning: pacote 'sp' foi compilado no R versão 4.4.3

## Warning: pacote 'Rcpp' foi compilado no R versão 4.4.3

shapefile_path_roubo <- "C:/Users/Vivian - H2R/Downloads/mba/git/mba_arrumado/nova abordagem/h3/h3_roubo_merge/content/h3_roubo_merge.shp"
shapefile_path_drogas <- "C:/Users/Vivian - H2R/Downloads/mba/git/mba_arrumado/nova abordagem/h3/h3_merge_drogas/content/h3_merge_drogas.shp" 

h3_roubo <- st_read(shapefile_path_roubo) %>%
  st_transform(crs = 31983)

## Reading layer `h3_roubo_merge' from data source 
##   `C:\Users\Vivian - H2R\Downloads\mba\git\mba_arrumado\nova abordagem\h3\h3_roubo_merge\content\h3_roubo_merge.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 297 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -46.7526 ymin: -23.67879 xmax: -46.57465 ymax: -23.50028
## Geodetic CRS:  WGS 84

h3_drogas <- st_read(shapefile_path_drogas) %>%
  st_transform(crs = 31983)

## Reading layer `h3_merge_drogas' from data source 
##   `C:\Users\Vivian - H2R\Downloads\mba\git\mba_arrumado\nova abordagem\h3\h3_merge_drogas\content\h3_merge_drogas.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 209 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -46.74932 ymin: -23.67879 xmax: -46.57389 ymax: -23.50449
## Geodetic CRS:  WGS 84

print("Colunas em h3_roubo:")

## [1] "Colunas em h3_roubo:"

print(colnames(h3_roubo))

## [1] "h3_index" "contagem" "geometry"

print("Colunas em h3_drogas:")

## [1] "Colunas em h3_drogas:"

print(colnames(h3_drogas))

## [1] "h3_index" "contagem" "geometry"

# Assegurar que ambos os data frames espaciais estão no mesmo CRS
h3_roubo <- st_transform(h3_roubo, crs = st_crs(h3_drogas))

# Encontrar os vizinhos mais próximos
nearest_idx <- st_nearest_feature(h3_roubo, h3_drogas)

# Unir os dados com base nos índices dos vizinhos mais próximos
h3_roubo_drogas <- cbind(h3_roubo, h3_drogas[nearest_idx,])

# Agora, remova a coluna de geometria duplicada (se necessário)
h3_roubo_drogas <- h3_roubo_drogas %>%
  select(-geometry.1)  # Remova a coluna de geometria duplicada

# Visualize as primeiras linhas para verificar
head(h3_roubo_drogas)

# *** Diagnóstico importante: Verificar NAs antes de converter para Spatial ***
print(paste("Valores NA em contagem (roubo) antes da conversão:", sum(is.na(h3_roubo_drogas$contagem)))) # contagem original de roubo

## [1] "Valores NA em contagem (roubo) antes da conversão: 0"

print(paste("Valores NA em contagem (drogas) antes da conversão:", sum(is.na(h3_roubo_drogas$contagem.1)))) # contagem da junção

## [1] "Valores NA em contagem (drogas) antes da conversão: 0"

# Conversão e limpeza ANTES de converter para Spatial
h3_roubo_drogas <- h3_roubo_drogas %>%
  mutate(contagem = as.numeric(as.character(contagem)),
         contagem.1 = as.numeric(as.character(contagem.1))) %>%
  filter(!is.na(contagem) & !is.na(contagem.1) & !is.nan(contagem) & !is.nan(contagem.1) & !is.infinite(contagem) & !is.infinite(contagem.1))

# Preparar os dados para GWR (CONVERTER PARA SPATIAL DEPOIS DE LIMPAR)
h3_roubo_drogas_sp <- as(h3_roubo_drogas, "Spatial")

# Verificações adicionais (agora DEVE funcionar)
print(class(h3_roubo_drogas_sp$contagem))

## [1] "numeric"

print(class(h3_roubo_drogas_sp$contagem.1))

## [1] "numeric"

print(paste("Número de linhas ANTES da filtragem:", nrow(h3_roubo_drogas_sp)))

## [1] "Número de linhas ANTES da filtragem: 297"

# Certifique-se de que 'contagem.x' e 'contagem.y' são variáveis numéricas (redundante, mas seguro)
h3_roubo_drogas_sp$contagem.x <- as.numeric(h3_roubo_drogas_sp$contagem)  # Use a coluna original "contagem"
h3_roubo_drogas_sp$contagem.y <- as.numeric(h3_roubo_drogas_sp$contagem.1) # Use a coluna junta "contagem.1"
# Verificar se há valores NA nas variáveis usadas na regressão
print(paste("Valores NA em contagem.x:", sum(is.na(h3_roubo_drogas_sp$contagem.x))))

## [1] "Valores NA em contagem.x: 0"

print(paste("Valores NaN em contagem.x:", sum(is.nan(h3_roubo_drogas_sp$contagem.x))))

## [1] "Valores NaN em contagem.x: 0"

print(paste("Valores Inf em contagem.x:", sum(is.infinite(h3_roubo_drogas_sp$contagem.x))))

## [1] "Valores Inf em contagem.x: 0"

print(paste("Valores NA em contagem.y:", sum(is.na(h3_roubo_drogas_sp$contagem.y))))

## [1] "Valores NA em contagem.y: 0"

print(paste("Valores NaN em contagem.y:", sum(is.nan(h3_roubo_drogas_sp$contagem.y))))

## [1] "Valores NaN em contagem.y: 0"

print(paste("Valores Inf em contagem.y:", sum(is.infinite(h3_roubo_drogas_sp$contagem.y))))

## [1] "Valores Inf em contagem.y: 0"

# Preparar os dados para GWR
h3_roubo_drogas_sp <- as(h3_roubo_drogas, "Spatial")

# Certifique-se de que 'contagem.x' e 'contagem.y' são variáveis numéricas
h3_roubo_drogas_sp$contagem.x <- as.numeric(h3_roubo_drogas_sp$contagem)  # Use a coluna original "contagem"
h3_roubo_drogas_sp$contagem.y <- as.numeric(h3_roubo_drogas_sp$contagem.1) # Use a coluna junta "contagem.1"

# Definir a largura de banda adaptativa usando cross-validation
# A largura de banda (bandwidth) é um parâmetro crucial na GWR. Ela determina o tamanho da janela espacial usada para ponderar as observações vizinhas.
# Uma largura de banda adaptativa permite que o tamanho da janela varie espacialmente, ajustando-se à densidade dos dados.
# Cross-validation (CV) é uma técnica para selecionar a largura de banda que minimiza o erro de previsão.
bw_gwr <- bw.gwr(contagem.x ~ contagem.y, 
                 data = h3_roubo_drogas_sp, 
                 approach = "CV", 
                 kernel = "gaussian", 
                 adaptive = TRUE)

## Adaptive bandwidth: 191 CV score: 20101202 
## Adaptive bandwidth: 126 CV score: 19272300 
## Adaptive bandwidth: 85 CV score: 18321155 
## Adaptive bandwidth: 60 CV score: 17429110 
## Adaptive bandwidth: 44 CV score: 16417585 
## Adaptive bandwidth: 35 CV score: 15469166 
## Adaptive bandwidth: 28 CV score: 14711061 
## Adaptive bandwidth: 25 CV score: 14449676 
## Adaptive bandwidth: 22 CV score: 14352628 
## Adaptive bandwidth: 21 CV score: 14266644 
## Adaptive bandwidth: 19 CV score: 12989039 
## Adaptive bandwidth: 19 CV score: 12989039

print(paste("Tamanho da banda utilizado:", bw_gwr))  # Imprime o tamanho da banda

## [1] "Tamanho da banda utilizado: 19"

# Executar o modelo GWR
# A função gwr.basic executa o modelo GWR.
# O parâmetro 'contagem.x ~ contagem.y' especifica a fórmula do modelo, onde 'contagem.x' é a variável dependente e 'contagem.y' é a variável independente.
# O parâmetro 'bw' especifica a largura de banda calculada usando bw.gwr.
# O parâmetro 'kernel' especifica a função kernel usada para ponderar as observações vizinhas. O kernel gaussiano é uma escolha comum.
# Os parâmetros 'hatmatrix = TRUE' e 'se.fit = TRUE' solicitam o cálculo da matriz hat e dos erros padrão dos coeficientes, respectivamente.
gwr_model <- gwr.basic(contagem.x ~ contagem.y, 
                       data = h3_roubo_drogas_sp, 
                       bw = bw_gwr, 
                       kernel = "gaussian", 
                       adaptive = TRUE)

# Imprimir os resultados do modelo GWR
# Imprimir os resultados do modelo GWR fornece informações sobre os coeficientes locais, os erros padrão, os valores t e os valores p.
# Essas informações podem ser usadas para avaliar a significância estatística da relação entre as variáveis e para identificar áreas onde a relação é mais forte ou mais fraca.
print(gwr_model)

##    ***********************************************************************
##    *                       Package   GWmodel                             *
##    ***********************************************************************
##    Program starts at: 2025-04-07 17:35:45.843401 
##    Call:
##    gwr.basic(formula = contagem.x ~ contagem.y, data = h3_roubo_drogas_sp, 
##     bw = bw_gwr, kernel = "gaussian", adaptive = TRUE)
## 
##    Dependent (y) variable:  contagem.x
##    Independent variables:  contagem.y
##    Number of data points: 297
##    ***********************************************************************
##    *                    Results of Global Regression                     *
##    ***********************************************************************
## 
##    Call:
##     lm(formula = formula, data = data)
## 
##    Residuals:
##     Min      1Q  Median      3Q     Max 
## -417.23 -128.51  -87.51   -0.17 1755.22 
## 
##    Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
##    (Intercept) 158.2856    15.9126   9.947  < 2e-16 ***
##    contagem.y    2.2203     0.5548   4.002 7.95e-05 ***
## 
##    ---Significance stars
##    Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
##    Residual standard error: 263.9 on 295 degrees of freedom
##    Multiple R-squared: 0.0515
##    Adjusted R-squared: 0.04828 
##    F-statistic: 16.02 on 1 and 295 DF,  p-value: 7.947e-05 
##    ***Extra Diagnostic information
##    Residual sum of squares: 20545243
##    Sigma(hat): 263.9032
##    AIC:  4158.739
##    AICc:  4158.821
##    BIC:  3889.901
##    ***********************************************************************
##    *          Results of Geographically Weighted Regression              *
##    ***********************************************************************
## 
##    *********************Model calibration information*********************
##    Kernel function: gaussian 
##    Adaptive bandwidth: 19 (number of nearest neighbours)
##    Regression points: the same locations as observations are used.
##    Distance metric: Euclidean distance metric is used.
## 
##    ****************Summary of GWR coefficient estimates:******************
##                    Min.   1st Qu.    Median   3rd Qu.   Max.
##    Intercept   40.58149  81.38859 125.71233 232.65804 573.90
##    contagem.y  -1.69096   0.36816   1.82600   3.69272  14.06
##    ************************Diagnostic information*************************
##    Number of data points: 297 
##    Effective number of parameters (2trace(S) - trace(S'S)): 20.95615 
##    Effective degrees of freedom (n-2trace(S) + trace(S'S)): 276.0438 
##    AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 4023.348 
##    AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 4004.681 
##    BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 3777.045 
##    Residual sum of squares: 11874133 
##    R-square value:  0.4518128 
##    Adjusted R-square value:  0.4100453 
## 
##    ***********************************************************************
##    Program stops at: 2025-04-07 17:35:45.864797

# Converter os resultados do GWR de volta para um objeto sf
# Os resultados do modelo GWR são armazenados em um objeto SDF (Spatial Data Frame).
# Para facilitar o mapeamento, é útil converter o SDF de volta para um objeto sf.
gwr_results_sf <- st_as_sf(gwr_model$SDF)

# Adicionar os coeficientes locais ao objeto sf
# Os coeficientes locais representam a relação entre as variáveis em cada localidade.
# Adicionar os coeficientes locais ao objeto sf permite mapeá-los e visualizar a variação espacial da relação.
h3_roubo_drogas$coef_drogas <- gwr_results_sf$contagem.y

# Criar mapas dos Coeficientes Locais (β), Valores Previstos (ŷ) e Resíduos
# tmap é um pacote R para criar mapas temáticos.
# tmap_mode("view") define o modo de visualização para interativo, permitindo explorar os mapas com zoom e pan.
tmap_mode("view")

## ℹ tmap mode set to "view".

# Extrair os valores previstos e resíduos do objeto gwr_model
h3_roubo_drogas$coef_drogas <- gwr_model$SDF$contagem.y #ja existia
h3_roubo_drogas$pred <- gwr_model$SDF$yhat  # Adicione os valores previstos
h3_roubo_drogas$residuals <- gwr_model$SDF$residual  # Adicione os resíduos

# Mapeamento dos resultados do GWR
tmap_mode("view")

## ℹ tmap mode set to "view".

# Mapa dos Coeficientes Locais
map_coef <- tm_shape(h3_roubo_drogas) +
  tm_fill("coef_drogas", title. = "Coeficiente Local de Drogas") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Coeficientes Locais de Drogas")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Mapa dos Valores Previstos
map_pred <- tm_shape(h3_roubo_drogas) +
  tm_fill("pred", title. = "Valores Previstos de Roubos") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Valores Previstos de Roubos")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Mapa dos Resíduos
map_resid <- tm_shape(h3_roubo_drogas) +
  tm_fill("residuals", title. = "Resíduos") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Resíduos")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Visualizar os mapas
map_coef

## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)

map_pred

map_resid

## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)

# Converter de volta para sf para usar st_write
h3_roubo_drogas_sf <- st_as_sf(h3_roubo_drogas)

# Salvar o shapefile com os resultados
st_write(h3_roubo_drogas_sf, "h3_roubo_drogas_gwr.shp", delete_layer = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## Deleting layer `h3_roubo_drogas_gwr' using driver `ESRI Shapefile'
## Writing layer `h3_roubo_drogas_gwr' to data source 
##   `h3_roubo_drogas_gwr.shp' using driver `ESRI Shapefile'
## Writing 297 features with 7 fields and geometry type Polygon.

# **Gráficos de Dispersão**

# 1. Valores Previstos vs. Observados
ggplot(h3_roubo_drogas, aes(x = pred, y = contagem)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # Linha de 45 graus
  labs(title = "Valores Previstos vs. Observados",
       x = "Valores Previstos (GWR)",
       y = "Valores Observados (Contagem de Roubos)") +
  theme_bw()

# 2. Resíduos vs. Valores Previstos
ggplot(h3_roubo_drogas, aes(x = pred, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") + # Linha horizontal em zero
  labs(title = "Resíduos vs. Valores Previstos",
       x = "Valores Previstos (GWR)",
       y = "Resíduos") +
  theme_bw()

# 3. Coeficientes Locais vs. Contagem de Drogas
ggplot(h3_roubo_drogas, aes(x = contagem.1, y = coef_drogas)) +
  geom_point() +
  labs(title = "Coeficientes Locais vs. Contagem de Drogas",
       x = "Contagem de Drogas",
       y = "Coeficientes Locais (Drogas)") +
  theme_bw()

# 4. Resíduos vs. Coordenadas Espaciais
# Calcular os centroides das geometrias
h3_roubo_drogas$centroid_x <- st_coordinates(st_centroid(h3_roubo_drogas$geometry))[,1]
h3_roubo_drogas$centroid_y <- st_coordinates(st_centroid(h3_roubo_drogas$geometry))[,2]

# Gráfico de Resíduos vs. Coordenada X
ggplot(h3_roubo_drogas, aes(x = centroid_x, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Resíduos vs. Coordenada X",
       x = "Coordenada X",
       y = "Resíduos") +
  theme_bw()

# Gráfico de Resíduos vs. Coordenada Y
ggplot(h3_roubo_drogas, aes(x = centroid_y, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Resíduos vs. Coordenada Y",
       x = "Coordenada Y",
       y = "Resíduos") +
  theme_bw()

# Este notebook R realiza uma análise de regressão ponderada geograficamente (GWR),
# cria gráficos de dispersão para avaliar os resultados do modelo e realiza uma regressão global
# para comparar os coeficientes estimados.

# Função para verificar e instalar pacotes
install_if_missing <- function(packages) {
  new_packages <- packages[!packages %in% installed.packages()[, "Package"]]
  if (length(new_packages) > 0) {
    install.packages(new_packages)
  }
}

# Lista de pacotes necessários
necessary_packages <- c("sf", "spdep", "tidyverse", "ncf", "tmap", "GWmodel", "dplyr", "ggplot2")
install_if_missing(necessary_packages)

# Carregar os pacotes
suppressPackageStartupMessages({
  library(sf)
  library(spdep)
  library(tidyverse)
  library(ncf)
  library(tmap)
  library(GWmodel)
  library(dplyr)
  library(ggplot2) # Garante que o ggplot2 está carregado
})

# Definir os caminhos para os shapefiles
shapefile_path_roubo <- "C:/Users/Vivian - H2R/Downloads/mba/git/mba_arrumado/nova abordagem/h3/h3_roubo_merge/content/h3_roubo_merge.shp"
shapefile_path_drogas <- "C:/Users/Vivian - H2R/Downloads/mba/git/mba_arrumado/nova abordagem/h3/h3_merge_drogas/content/h3_merge_drogas.shp"

# Ler os shapefiles
h3_roubo <- st_read(shapefile_path_roubo) %>%
  st_transform(crs = 31983)

## Reading layer `h3_roubo_merge' from data source 
##   `C:\Users\Vivian - H2R\Downloads\mba\git\mba_arrumado\nova abordagem\h3\h3_roubo_merge\content\h3_roubo_merge.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 297 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -46.7526 ymin: -23.67879 xmax: -46.57465 ymax: -23.50028
## Geodetic CRS:  WGS 84

h3_drogas <- st_read(shapefile_path_drogas) %>%
  st_transform(crs = 31983)

## Reading layer `h3_merge_drogas' from data source 
##   `C:\Users\Vivian - H2R\Downloads\mba\git\mba_arrumado\nova abordagem\h3\h3_merge_drogas\content\h3_merge_drogas.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 209 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -46.74932 ymin: -23.67879 xmax: -46.57389 ymax: -23.50449
## Geodetic CRS:  WGS 84

# Assegurar que ambos os data frames espaciais estão no mesmo CRS
h3_roubo <- st_transform(h3_roubo, crs = st_crs(h3_drogas))

# Encontrar os vizinhos mais próximos
nearest_idx <- st_nearest_feature(h3_roubo, h3_drogas)

# Unir os dados com base nos índices dos vizinhos mais próximos
h3_roubo_drogas <- cbind(h3_roubo, h3_drogas[nearest_idx,])

# Agora, remova a coluna de geometria duplicada (se necessário)
h3_roubo_drogas <- h3_roubo_drogas %>%
  select(-geometry.1)  # Remova a coluna de geometria duplicada

# Conversão e limpeza ANTES de converter para Spatial
h3_roubo_drogas <- h3_roubo_drogas %>%
  mutate(contagem = as.numeric(as.character(contagem)),
         contagem.1 = as.numeric(as.character(contagem.1))) %>%
  filter(!is.na(contagem) & !is.na(contagem.1) & !is.nan(contagem) & !is.nan(contagem.1) & !is.infinite(contagem) & !is.infinite(contagem.1))

# Preparar os dados para GWR (CONVERTER PARA SPATIAL DEPOIS DE LIMPAR)
h3_roubo_drogas_sp <- as(h3_roubo_drogas, "Spatial")

# Certifique-se de que 'contagem.x' e 'contagem.y' são variáveis numéricas
h3_roubo_drogas_sp$contagem.x <- as.numeric(h3_roubo_drogas_sp$contagem)  # Use a coluna original "contagem"
h3_roubo_drogas_sp$contagem.y <- as.numeric(h3_roubo_drogas_sp$contagem.1) # Use a coluna junta "contagem.1"

# Definir a largura de banda adaptativa usando cross-validation
bw_gwr <- bw.gwr(contagem.x ~ contagem.y,
                 data = h3_roubo_drogas_sp,
                 approach = "CV",
                 kernel = "gaussian",
                 adaptive = TRUE)

## Adaptive bandwidth: 191 CV score: 20101202 
## Adaptive bandwidth: 126 CV score: 19272300 
## Adaptive bandwidth: 85 CV score: 18321155 
## Adaptive bandwidth: 60 CV score: 17429110 
## Adaptive bandwidth: 44 CV score: 16417585 
## Adaptive bandwidth: 35 CV score: 15469166 
## Adaptive bandwidth: 28 CV score: 14711061 
## Adaptive bandwidth: 25 CV score: 14449676 
## Adaptive bandwidth: 22 CV score: 14352628 
## Adaptive bandwidth: 21 CV score: 14266644 
## Adaptive bandwidth: 19 CV score: 12989039 
## Adaptive bandwidth: 19 CV score: 12989039

# Executar o modelo GWR
gwr_model <- gwr.basic(contagem.x ~ contagem.y,
                       data = h3_roubo_drogas_sp,
                       bw = bw_gwr,
                       kernel = "gaussian",
                       adaptive = TRUE)

# Extrair os valores previstos e resíduos do objeto gwr_model
h3_roubo_drogas$coef_drogas <- gwr_model$SDF$contagem.y  # já existia
h3_roubo_drogas$pred <- gwr_model$SDF$yhat  # Adicione os valores previstos
h3_roubo_drogas$residuals <- gwr_model$SDF$residual  # Adicione os resíduos

# **Regressão Global**
# A regressão global serve como um ponto de referência para comparar com os resultados do GWR.
# Ela assume que a relação entre roubos e contagem de drogas é constante em toda a área de estudo.
global_model <- lm(contagem ~ contagem.1, data = h3_roubo_drogas)
summary(global_model)

## 
## Call:
## lm(formula = contagem ~ contagem.1, data = h3_roubo_drogas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -417.23 -128.51  -87.51   -0.17 1755.22 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 158.2856    15.9126   9.947  < 2e-16 ***
## contagem.1    2.2203     0.5548   4.002 7.95e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 263.9 on 295 degrees of freedom
## Multiple R-squared:  0.0515, Adjusted R-squared:  0.04828 
## F-statistic: 16.02 on 1 and 295 DF,  p-value: 7.947e-05

# Extrair coeficientes da regressão global
global_coef <- coef(global_model)[2] # Coeficiente da contagem de drogas
global_intercept <- coef(global_model)[1] # Intercepto

# Adicionar valores ajustados da regressão global ao dataframe
h3_roubo_drogas$global_pred <- predict(global_model)

# Converter de volta para sf para usar st_write
h3_roubo_drogas_sf <- st_as_sf(h3_roubo_drogas)

# Salvar o shapefile com os resultados
st_write(h3_roubo_drogas_sf, "h3_roubo_drogas_gwr.shp", delete_layer = TRUE)

## Warning in abbreviate_shapefile_names(obj): Field names abbreviated for ESRI
## Shapefile driver

## Deleting layer `h3_roubo_drogas_gwr' using driver `ESRI Shapefile'
## Writing layer `h3_roubo_drogas_gwr' to data source 
##   `h3_roubo_drogas_gwr.shp' using driver `ESRI Shapefile'
## Writing 297 features with 8 fields and geometry type Polygon.

# **Gráficos de Dispersão**

# 1. Valores Previstos (GWR) vs. Observados
# Este gráfico ajuda a avaliar o quão bem o modelo GWR está prevendo os valores reais de roubo.
# Idealmente, os pontos devem se agrupar próximos à linha de 45 graus.
ggplot(h3_roubo_drogas, aes(x = pred, y = contagem)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # Linha de 45 graus
  labs(title = "Valores Previstos (GWR) vs. Observados",
       x = "Valores Previstos (GWR)",
       y = "Valores Observados (Contagem de Roubos)") +
  theme_bw()

# 2. Valores Previstos (Global) vs. Observados
# Similar ao gráfico anterior, mas para o modelo de regressão global.
# Comparar este gráfico com o anterior ajuda a determinar se o GWR oferece uma melhoria em relação ao modelo global.
ggplot(h3_roubo_drogas, aes(x = global_pred, y = contagem)) +
  geom_point() +
  geom_abline(intercept = 0, slope = 1, color = "red", linetype = "dashed") + # Linha de 45 graus
  labs(title = "Valores Previstos (Global) vs. Observados",
       x = "Valores Previstos (Global)",
       y = "Valores Observados (Contagem de Roubos)") +
  theme_bw()

# 3. Coeficientes Locais vs. Contagem de Drogas
# Este gráfico examina a relação entre os coeficientes locais estimados pelo GWR e a contagem de drogas.
# Ele pode revelar padrões interessantes sobre como a influência da contagem de drogas nos roubos varia espacialmente.
# Por exemplo, pode haver uma relação não linear ou uma tendência de aumento/diminuição do coeficiente com o aumento da contagem de drogas.
ggplot(h3_roubo_drogas, aes(x = contagem.1, y = coef_drogas)) +
  geom_point() +
  labs(title = "Coeficientes Locais (GWR) vs. Contagem de Drogas",
       x = "Contagem de Drogas",
       y = "Coeficientes Locais (GWR)") +
  theme_bw()

# 4. Resíduos (GWR) vs. Valores Previstos (GWR)
# Este gráfico ajuda a verificar a homocedasticidade (variância constante dos erros) e a linearidade do modelo GWR.
# Idealmente, os resíduos devem ser distribuídos aleatoriamente em torno de zero, sem padrões óbvios.
# Padrões como um funil ou uma curva podem indicar problemas com o modelo.
ggplot(h3_roubo_drogas, aes(x = pred, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") + # Linha horizontal em zero
  labs(title = "Resíduos (GWR) vs. Valores Previstos (GWR)",
       x = "Valores Previstos (GWR)",
       y = "Resíduos (GWR)") +
  theme_bw()

# 5. Resíduos vs. Coordenadas Espaciais
# Calcular os centroides das geometrias
h3_roubo_drogas$centroid_x <- st_coordinates(st_centroid(h3_roubo_drogas$geometry))[,1]
h3_roubo_drogas$centroid_y <- st_coordinates(st_centroid(h3_roubo_drogas$geometry))[,2]

# Gráfico de Resíduos vs. Coordenada X
# Estes gráficos ajudam a verificar se há padrões espaciais nos resíduos.
# Se houver um padrão, isso pode indicar que o modelo não está capturando toda a variação espacial nos dados.
ggplot(h3_roubo_drogas, aes(x = centroid_x, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Resíduos (GWR) vs. Coordenada X",
       x = "Coordenada X",
       y = "Resíduos (GWR)") +
  theme_bw()

# Gráfico de Resíduos vs. Coordenada Y
ggplot(h3_roubo_drogas, aes(x = centroid_y, y = residuals)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  labs(title = "Resíduos (GWR) vs. Coordenada Y",
       x = "Coordenada Y",
       y = "Resíduos (GWR)") +
  theme_bw()

# Mapeamento dos resultados do GWR
tmap_mode("view")

## ℹ tmap mode set to "view".

# Mapa dos Coeficientes Locais
# Este mapa mostra como o coeficiente da contagem de drogas varia espacialmente.
# Áreas com coeficientes mais altos indicam onde a contagem de drogas tem uma influência maior nos roubos.
map_coef <- tm_shape(h3_roubo_drogas) +
  tm_fill("coef_drogas", title. = "Coeficiente Local de Drogas") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Coeficientes Locais de Drogas")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Mapa dos Valores Previstos (GWR)
# Este mapa mostra os valores previstos de roubos pelo modelo GWR.
# Ele pode ser usado para identificar áreas com alto risco de roubo.
map_pred <- tm_shape(h3_roubo_drogas) +
  tm_fill("pred", title. = "Valores Previstos de Roubos (GWR)") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Valores Previstos de Roubos (GWR)")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Mapa dos Resíduos (GWR)
# Este mapa mostra os resíduos do modelo GWR.
# Ele pode ser usado para identificar áreas onde o modelo está subestimando ou superestimando os roubos.
map_resid <- tm_shape(h3_roubo_drogas) +
  tm_fill("residuals", title. = "Resíduos (GWR)") + # correção para tmap v4
  tm_borders() +
  tm_layout(title = "Mapa dos Resíduos (GWR)")

## [tm_fill()] Argument `title.` unknown.
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`

# Visualizar os mapas
map_coef

## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)

map_pred

map_resid

## [scale] tm_polygons:() the data variable assigned to 'fill' contains positive and negative values, so midpoint is set to 0. Set 'midpoint = NA' in 'fill.scale = tm_scale_intervals(<HERE>)' to use all visual values (e.g. colors)